#!/usr/bin/env julia
"""
Adaptive Refinement Demonstration

This example demonstrates the adaptive p-refinement capabilities of NSEMSolver.jl
for L-shaped domain problems with corner singularities.
"""

using NSEMSolver

function demonstrate_lshape_refinement()
    println("🔄 L-Shaped Domain with Adaptive Refinement")
    println("=" * 50)
    
    # Create options with adaptive refinement enabled
    options = NSOptions(
        N = 3,                    # Start with lower polynomial order
        n_block = 4,             # More blocks for L-shaped domain
        nu = 0.005,              # Lower viscosity for interesting flow
        tfinal = 2.0,            # Longer time to see adaptation
        cfl = 0.4,
        tol = 0.01,              # Relaxed tolerance for demonstration
        solver = :julia,
        method = :sem,
        adaptive_refinement = true,  # Enable adaptive refinement
        refinement_levels = 4,   # Multiple refinement levels
        verbose = true,
        save_history = true
    )
    
    println("🏗️  Domain Configuration:")
    println("  L-shaped domain with corner singularity")
    println("  Adaptive p-refinement: enabled")
    println("  Refinement levels: $(options.refinement_levels)")
    println("  Base polynomial order: $(options.N)")
    println()
    
    # Solve with adaptive refinement
    println("🚀 Starting adaptive refinement simulation...")
    result = solve_navier_stokes_2d(options)
    
    # Display results
    println("\n📊 Adaptive Refinement Results:")
    println("  Converged: $(result.converged)")
    println("  Total iterations: $(result.iterations)")
    println("  Final residual: $(@sprintf("%.2e", result.residual_norm))")
    println("  Solve time: $(@sprintf("%.3f", result.solve_time)) seconds")
    
    # Analyze refinement map
    if result.multidomain !== nothing && result.multidomain.refinement_map !== nothing
        println("\n🔍 Refinement Analysis:")
        refinement_map = result.multidomain.refinement_map
        matrix_dist = result.multidomain.matrix_dist
        
        # Count elements at each refinement level
        refinement_counts = Dict{Int, Int}()
        total_active = 0
        
        for idx in CartesianIndices(refinement_map)
            if matrix_dist[idx] == 1  # Active domain element
                level = refinement_map[idx]
                refinement_counts[level] = get(refinement_counts, level, 0) + 1
                total_active += 1
            end
        end
        
        println("  Total active elements: $total_active")
        for level in sort(collect(keys(refinement_counts)))
            count = refinement_counts[level]
            percentage = (count / total_active) * 100
            p_order = get_polynomial_order_from_refinement(level)
            println("  Level $level (p=$p_order): $count elements ($(@sprintf("%.1f", percentage))%)")
        end
        
        # Find maximum refinement locations
        max_level = maximum(refinement_map[matrix_dist .== 1])
        max_indices = findall(x -> x == max_level, refinement_map)
        println("  Maximum refinement level: $max_level")
        println("  High refinement locations: $(length(max_indices)) elements")
        
        # Show corner refinement effectiveness
        corner_points = result.multidomain.corner_points
        if length(corner_points) > 0
            println("  Corner point refinement:")
            for (i, corner) in enumerate(corner_points)
                if length(corner) == 2
                    ci, cj = corner
                    level = refinement_map[ci, cj]
                    p_order = get_polynomial_order_from_refinement(level)
                    println("    Corner $i at ($ci, $cj): level $level (p=$p_order)")
                end
            end
        end
    end
    
    return result
end

function compare_adaptive_vs_uniform()
    println("\n⚖️  Adaptive vs Uniform Refinement Comparison")
    println("-" * 50)
    
    base_options = NSOptions(
        N = 3,
        n_block = 3,
        nu = 0.01,
        tfinal = 1.0,
        cfl = 0.4,
        tol = 1e-5,
        solver = :julia,
        verbose = false
    )
    
    # Run with uniform refinement (disabled adaptive)
    println("🔧 Running uniform refinement case...")
    uniform_options = NSOptions(base_options; adaptive_refinement = false)
    uniform_result = solve_navier_stokes_2d(uniform_options)
    
    # Run with adaptive refinement
    println("🔄 Running adaptive refinement case...")
    adaptive_options = NSOptions(base_options; adaptive_refinement = true, refinement_levels = 3)
    adaptive_result = solve_navier_stokes_2d(adaptive_options)
    
    # Compare results
    println("\n📊 Comparison Results:")
    println("  Uniform Refinement:")
    println("    Converged: $(uniform_result.converged)")
    println("    Iterations: $(uniform_result.iterations)")
    println("    Time: $(@sprintf("%.3f", uniform_result.solve_time)) s")
    println("    Final residual: $(@sprintf("%.2e", uniform_result.residual_norm))")
    
    println("  Adaptive Refinement:")
    println("    Converged: $(adaptive_result.converged)")
    println("    Iterations: $(adaptive_result.iterations)")
    println("    Time: $(@sprintf("%.3f", adaptive_result.solve_time)) s")
    println("    Final residual: $(@sprintf("%.2e", adaptive_result.residual_norm))")
    
    # Calculate improvement metrics
    if uniform_result.converged && adaptive_result.converged
        time_ratio = uniform_result.solve_time / adaptive_result.solve_time
        iter_ratio = uniform_result.iterations / adaptive_result.iterations
        
        println("\n🏆 Performance Improvement:")
        if adaptive_result.solve_time < uniform_result.solve_time
            println("  ⚡ Adaptive is $(@sprintf("%.1fx", time_ratio)) faster")
        else
            println("  ⏳ Uniform is $(@sprintf("%.1fx", 1/time_ratio)) faster")
        end
        
        if adaptive_result.iterations < uniform_result.iterations
            println("  🎯 Adaptive uses $(@sprintf("%.1fx", iter_ratio)) fewer iterations")
        else
            println("  🔄 Uniform uses $(@sprintf("%.1fx", 1/iter_ratio)) fewer iterations")
        end
        
        # Accuracy comparison
        if adaptive_result.residual_norm < uniform_result.residual_norm
            residual_improvement = uniform_result.residual_norm / adaptive_result.residual_norm
            println("  🎯 Adaptive achieves $(@sprintf("%.1fx", residual_improvement)) better accuracy")
        end
    end
    
    return uniform_result, adaptive_result
end

function demonstrate_error_indicators()
    println("\n🧠 Physics-Based Error Indicators")
    println("-" * 40)
    
    # Create a simple test case
    options = NSOptions(
        N = 4,
        n_block = 3,
        nu = 0.01,
        tfinal = 0.5,
        tol = 1e-6,
        adaptive_refinement = true,
        verbose = false
    )
    
    println("🔬 Testing error indicator calculation...")
    
    result = solve_navier_stokes_2d(options)
    
    if result.converged && result.multidomain !== nothing
        # Extract velocity fields for error analysis
        u = result.u
        v = result.v
        
        # Create dummy SEM operators for error calculation
        sem_operators = create_sem_operators(options.N, options.n_block, 1.0, 2)
        D1, D2 = sem_operators.D1, sem_operators.D2
        
        # Calculate error indicators
        error_indicators = calculate_physics_error_indicators(
            u, v, nothing, D1, D2, options.nu, 2
        )
        
        # Analyze error distribution
        error_stats = analyze_error_distribution(error_indicators, result.multidomain.matrix_dist)
        
        println("📈 Error Indicator Statistics:")
        println("  Maximum error: $(@sprintf("%.3e", error_stats.max_error))")
        println("  Mean error: $(@sprintf("%.3e", error_stats.mean_error))")
        println("  Standard deviation: $(@sprintf("%.3e", error_stats.std_error))")
        println("  High error regions: $(error_stats.high_error_count) elements")
        
        # Show spatial distribution
        high_error_threshold = error_stats.mean_error + 2 * error_stats.std_error
        println("  High error threshold: $(@sprintf("%.3e", high_error_threshold))")
        
        return error_indicators, error_stats
    else
        println("❌ Could not compute error indicators (simulation failed)")
        return nothing, nothing
    end
end

function analyze_error_distribution(error_indicators::Array, matrix_dist::Array)
    # Extract active domain errors
    active_errors = error_indicators[matrix_dist .== 1]
    
    max_error = maximum(active_errors)
    mean_error = sum(active_errors) / length(active_errors)
    std_error = sqrt(sum((active_errors .- mean_error).^2) / (length(active_errors) - 1))
    
    # Count high error regions
    threshold = mean_error + std_error
    high_error_count = count(x -> x > threshold, active_errors)
    
    return (
        max_error = max_error,
        mean_error = mean_error,
        std_error = std_error,
        high_error_count = high_error_count,
        total_elements = length(active_errors)
    )
end

function main()
    println("🔄 NSEMSolver.jl - Adaptive Refinement Demonstration")
    println("=" * 60)
    
    # Demo 1: L-shaped domain with adaptive refinement
    lshape_result = demonstrate_lshape_refinement()
    
    # Demo 2: Compare adaptive vs uniform refinement
    uniform_result, adaptive_result = compare_adaptive_vs_uniform()
    
    # Demo 3: Demonstrate error indicators
    error_indicators, error_stats = demonstrate_error_indicators()
    
    # Summary
    println("\n" * "=" * 60)
    println("🎯 Adaptive Refinement Summary:")
    
    if lshape_result.converged
        println("✅ L-shaped domain simulation: SUCCESS")
    else
        println("❌ L-shaped domain simulation: FAILED")
    end
    
    if adaptive_result.converged && uniform_result.converged
        time_improvement = uniform_result.solve_time / adaptive_result.solve_time
        if time_improvement > 1.0
            println("⚡ Adaptive refinement achieved $(@sprintf("%.1fx", time_improvement)) speedup")
        else
            println("⏳ Uniform refinement was $(@sprintf("%.1fx", 1/time_improvement)) faster")
        end
    end
    
    if error_stats !== nothing
        println("🧠 Error indicators successfully computed")
        println("📊 High error regions identified: $(error_stats.high_error_count) elements")
    end
    
    println("\n🎉 Adaptive refinement demonstration completed!")
    
    return lshape_result, adaptive_result, uniform_result
end

# Run the demonstration
if abspath(PROGRAM_FILE) == @__FILE__
    main()
end